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Oh Abstract 
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In this paper we use the Generalized Multiscale Finite Element Method (GMsFEM) framework, 

introduced in [|20ll , in order to solve nonlinear elliptic equations with high-contrast coefficients. 

The proposed solution method involves linearizing the equation so that coarse-grid quantities of 

previous solution iterates can be regarded as auxiliary parameters within the problem formulation. 

> With this convention, we systematically construct respective coarse solution spaces that lend them- 

00 

selves to either continuous Galerkin (CG) or discontinuous Galerkin (DG) global formulations. 
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Here, we use Symmetric Interior Penalty Discontinuous Galerkin approach. Both methods yield a 
predictable error decline that depends on the respective coarse space dimension, and we illustrate 
the effectiveness of the CG and DG formulations by offering a variety of numerical examples. 
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1. Introduction 
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Nonlinear partial differential equations represent a class of problems that have applications in 
many scientific communities lfT6ll35l . Forchheimer flow and nonlinear elasticity are two particular 
examples of physical processes that are modeled by nonlinear equations [fT4ll27ll . In addition to 
difficulties associated with the nonlinearity, these types of problems often involve coefficients that 
exhibit high-contrast, heterogeneous behavior. For example, when modeling subsurface flow, the 
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underlying permeability field is often represented by a high-contrast coefficient in the pressure 
equation. One approach for solving a high-contrast, nonlinear equation is to linearize the problem 
and use an iterative method for obtaining the solution. For example, a Picard iteration yields an 
iterative process where a previous solution iterate is directly used in order to update the solution at 
the current iteration. In this case, a final solution is obtained when a suitable tolerance between the 
current and previous iteration is reached. While relatively easy to implement, iterative techniques 
typically require a repeated number of solves in order to obtain a convergent solution. In the case 
of a nonlinear elliptic equation, each iteration requires the numerical solution of a large system of 
equations that depends on the previous iterate. Thus, computing solutions on a fully resolved mesh 
quickly becomes a prohibitively expensive task. As such, techniques that allow for a more efficient 
computational procedure with a suitable level of accuracy are desirable. 

The past few decades have seen the development of various multiscale solution techniques for 
capturing small scale effects on a coarse grid [HI El |24l I28l - l30ll . The multiscale finite element 
methods (MsFEM's) that we consider in this paper hinge on the construction of coarse spaces that 
are spanned by a set of independently computed multiscale basis functions. The multiscale basis 
functions are then coupled via a respective global formulation in order to compute the solution. 
In particular, solutions may be computed on a coarse grid while maintaining the fine-scale effects 
that are embedded into the basis functions. While standard multiscale methods have proven effec- 
tive for a variety of applications (see, e.g., Il23] - l25l l30l0 . in this paper we consider a more recent 
framework in which the coarse spaces may be systematically enriched to converge to the fine grid 
solution J51|2Tl|22l[34l|. More specifically, additional basis functions are chosen based on local- 
ized eigenvalue problems that capture the underlying behavior of the system. In this case, we may 
carefully choose the number of basis functions (and dimension of the coarse space) such that we 
achieve a desired level of accuracy. In this paper we additionally show that the systematic enrich- 
ment of coarse spaces is flexible with respect to the global formulation that is chosen to couple the 
resulting basis functions. 

To treat the nonlinear elliptic equation considered in this paper we make use of the Generalized 
Multiscale Finite Element Method (GMsFEM) which was introduced in ll20l . In order to do so, we 
apply a Picard iteration and treat an upscaled quantity of a previous solution iterate as a parameter 
in the problem. With this convention we follow an offline-online procedure in which the coarse 
space construction is split into two distinct stages; offline and online (see [[H OH [341 [371). The 
main goal of this approach is to allow for the efficient construction of an online space (and an 
online solution) for each fixed parameter value and iteration. In the process, we precompute a 
larger-dimensional, parameter-independent offline space that accounts for an appropriate range 



of parameter values that may be used in the online stage. As construction of the offline space 
will constitute a one-time preprocessing step, only the online space will require additional work 
within the solution procedure. In the offline stage we first choose a fixed set of parameter values 
and generate an associated set of "snapshot" functions by solving localized problems on specified 
coarse subdomains. The functions obtained through this step constitute a snapshot space which 
will be used in the offline space construction. To construct the offline space we solve localized 
eigenvalue problems that use averaged quantities of the parameter(s) of interest within the space 
of snapshots. We then keep a certain number of eigenfunctions (based on some criterion) to form 
the offline space. At the online stage we solve similar localized problems using a fixed parameter 
value within the offline space, and keep a certain number of eigenfunctions for the online space 
construction. 

In this paper we consider the continuous Galerkin (CG) and discontinuous Galerkin (DG) for- 
mulations for the global coupling of the online basis functions. We show that each method offers a 
suitable solution technique, however, at this point we highlight some distinguishing characteristics 
of the repsective methods as motivation for considering both formulations. For the nonlinear ellip- 
tic equation considered in this paper, the CG coupling yields a bilinear form that closely resembles 
the standard finite element method (FEM). In particular, the integrations that define the CG for- 
mulation are taken over the whole domain, and result in a reduced-order system of equations that 
is similar in nature to the fine-scale system. As such, the ease of implementation, classical FEM 
analogues, and well understood structure make CG a tractable method for coupling the coarse basis 
functions in order to solve the global problem [|28l . While the discontinuous Galerkin formulation 
is arguably more delicate than its CG counterpart, DG offers an attractive feature such as it does 
not require partition of unity functions to couple basis functions. Both methods are shown to be 
suitable coupling mechanisms within the GMsFEM framework that is described in this paper. In 
particular, an increase in the size of the online coarse space yields a predictable error decline, and 
the error is shown to behave according to previous error estimates that depend on the eigenvalue 
behavior. The flexibility of the coarse space enrichment, along with the choice of using CG or DG 
as the global coupling mechanism, makes GMsFEM a robust and suitable technique for solving 
the model equation that we consider in this paper. A variety of numerical examples are presented 
to validate the performance of the proposed method. 

We note that some numerical results for GMsFEM in the context of continuous Galerkin meth- 
ods for nonlinear equations are presented in ll20ll . These numerical results are mostly presented 
to demonstrate the main concepts of GMsFEM and we do not have careful studies for nonlinear 
problems in [20]. Moreover, the numerical results presented in [20] use reduced basis approach 



to identify dominant eigenmodes which is different from the local mode decomposition approach 
presented here. Moreover, the current paper also studies DG approach for nonlinear equations. 

The organization of the paper is as follows. In Sect. [2] we introduce the model problem, the 
iterative procedure, and notation to be used throughout the paper. In Sect. [3] we carefully describe 
the coarse space enrichment procedure, and introduce the continuous and discontinuous Galerkin 



global coupling formulations. In particular, Subsect. 3.1 is devoted to the offline-online coarse 



space construction, and in Subsect. 3.2 we describe the CG and DG global coupling procedures. A 
variety of numerical examples are presented in Sect. |4]to validate the performance of the proposed 
approaches, and in Sect. [5] we offer some concluding remarks. 

2. Preliminaries 

In this paper we consider non-linear, elliptic equations of the form 

-div(«(z;u)Vu) = /inD, (1) 

where u = on dD. We assume that u is bounded above and below, i.e., u < u(x) < un, where 
m and un are pre-defined constants. We will also assume that the interval [u , u^] is divided into 
iV equal regions whose endpoints are given by u < U\ < ... < Wat-i < %■ 
In order to solve Eq. ([I]) we will consider a Picard iteration 

- div(«(a;; u n (x)) Vu n+1 (x)) = f in D, (2) 

where superscripts involving n denote respective iteration levels. To discretize ([2]), we next in- 
troduce the notion of fine and coarse grids. We let T H be a usual conforming partition of the 
computational domain D into finite elements (triangles, quadrilaterals, tetrahedrals, etc.). We re- 
fer to this partition as the coarse grid and assume that each coarse subregion is partitioned into a 
connected union of fine grid blocks. The fine grid partition will be denoted by T h . We use {xj}^ 
(where N v the number of coarse nodes) to denote the vertices of the coarse mesh T H , and define 
the neighborhood of the node Xi by 



UJ; 



[JiKjET 11 ; XiEKj}. (3) 



See Fig.[T]for an illustration of neighborhoods and elements subordinated to the coarse discretiza- 
tion. We emphasize the use of U{ to denote a coarse neighborhood, and K to denote a coarse 
element throughout the paper. 
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Figure 1 : Illustration of a coarse neighborhood and coarse element 



Next, we briefly outline the global coupling and the role of coarse basis functions for the re- 
spective formulations that we consider. For the discontinuous Galerkin (DG) formulation, we 
will use a coarse element K as the support for basis functions, and for the continuous Galerkin 
(CG) formulation, we will use w* as the support of basis functions. To further motivate the coarse 
basis construction, we offer a brief outline of the global coupling associated with the CG formu- 
lation below. For the purposes of this description, we formally denote the CG basis functions 
by ip^. In particular, we note that the proposed approach will employ the use of multiple basis 
functions per coarse neighborhood. In turn, the CG solution at n-th iteration will be sought as 
u^(x; /i) = Y^,i k c k' l l J k 1 ( x '-> A 4 )' where ^(x; /i) are the basis functions for n-th iteration, and /i 
is used to denote dependence on the previous solution. We note that a main consideration of our 
method is to allow for rapid calculations of basis functions at each iteration. 

Once the basis functions are identified, the CG global coupling is given through the variational 
form 

a(u™, v- fi) = (f, v), for all v G V™, (4) 

where V^ G is used to denote the space formed by those basis functions. 

We also note that an appropriate set of basis functions defined on each coarse element K may 
be respectively coupled via a discontinuous Galerkin formulation (see e.g., [|4l [T5ll36l '). 



3. CG and DG GMsFEM for nonlinear problems 

3.1. Local basis functions 

To motivate the local basis construction, we first introduce an approximation to the solution of 
Eq. (Q given by 

- div(K(x;u n (x))Vu n+1 (x)) = f inD, (5) 

where u denotes the average of u in each coarse region (either K or Wj, depending on the desired 
formulation). Because the variation in u" is not known a priori, we use /i to represent the depen- 
dence of the solution on u n . As part of the iterative solution process, multiscale basis functions 
will be computed for a selected number of the parameter values at the offline stage, and we will 
compute multiscale basis functions for each new value of u n at the online stage. In this section we 
will describe these details, and note that we maintain the convention of denoting u by the param- 
eter [i. We omit the iterative index n (and n + 1) for additional notational brevity, although note 
that the iterative process should be clearly implied. 

With the notational conventions in place we now describe the offline- online computational 
procedure, and elaborate on some applicable choices for the associated bilinear forms to be used 
in the coarse space construction. Below we offer a general outline for the procedure. 

1. Offline computations: 

- 1 .0. Coarse grid generation. 

- 1.1. Construction of snapshot space that will be used to compute an offline space. 

- 1.2. Construction of a small dimensional offline space by performing dimension reduc- 
tion in the space of global snapshots. 

2. Online computations: 

- 2.1. For each input parameter, compute multiscale basis functions. 

- 2.2. Solution of a coarse-grid problem for any force term and boundary condition. 

- 2.3. Iterative solvers, if needed. 

In the offline computation, we first construct a snapshot space V^ ap , where r denotes either a 
coarse neighborhood cu,i in the continuous Galerkin case, or a coarse element K in the discontin- 
uous Galerkin case (refer back to Fig.[T|). Construction of the snapshot space involves solving the 
local problems for various choices of input parameters, and we describe the details below. 



In order to construct the space of snapshots we propose to solve the following eigenvalue 
problem on a coarse domain r: 

M^wr = xrsfaWiT in *> ^ 

where /jlj (j = 1, . . . , J) is a specified set of fixed parameter values, and we again emphasize that 
r denotes a different coarse subdomain (either a coarse neighborhood U{ or coarse element K) 
depending on whether we consider the CG or DG problem formulation. We are careful to note that 
zero Neumann boundary conditions are generally used to solve eigenvalue problem, except in the 
DG case when Dirichlet conditions are used on element boundaries that coincide with the global 
domain. The matrices in Eq. ([6]) are defined as 

A(hj) = [ain^mn] = / /t(x;/i.,)V0 n -V0 m and S(/j,j) = [s(/^-) mn ] = / k(x; fij)(j) n (j) m , O) 

where <\> n denotes the standard bilinear, fine-scale basis functions and k will be carefully introduced 
in the next section. We note that Eq. (|6]) is the discretized form of the continuous equation 

j- / / \t-t / T,snap\ \ r,snap / r.snap 

-div(«(x, fij) Vi/>£ v ) = \j v *Pu m t. 

For brevity of notation we now omit the superscript r for eigenvalue problems, yet it is assumed 
throughout this section that the offline and online space computations are localized to respective 
coarse subdomains. After solving Eq. (|6]), we keep the first L, eigenfunctions corresponding to the 
dominant eigenvalues (asymptotically vanishing in this case) to form the space 

^snap = span{^ s " ap : 1 < j < J and 1 < / < L { }, 

for each coarse neighborhood tOi (or coarse element K). 

We reorder the snapshot functions using a single index to create the matrix 



-*l-snap 



,snap ,snap 

V\ ' • • • ' m na „ 



where M snap denotes the total number of functions to keep in the snapshot matrix construction. 

In order to construct the offline space V^, we perform a dimension reduction of the space of 
snapshots using an auxiliary spectral decomposition. The main objective is to use the offline space 
to efficiently (and accurately) construct a set of multiscale basis functions for each // value in the 
online stage. More precisely, we seek a subspace of the snapshot space such that it can approximate 



any element of the snapshot space in the appropriate sense defined via auxiliary bilinear forms. At 
the offline stage the bilinear forms are chosen to be parameter-independent, such that there is no 
need to reconstruct the offline space for each \i value. The analysis in [22] motivates the following 
eigenvalue problem in the space of snapshots: 

A off Vf = \fS oS Vf, (8) 



where 



A°* = [O* ] = |^,/i)VC P " Wr P = ^ap^snap 



and 

*-> = l S mn\ = / / H X ;/ i ) ? /m Vn = -^snap^-Rsnap, 

where k(x, /Z), and k(x, JJl) are domain-based averaged coefficients with JL chosen as the average of 
pre-selected /Vs. We note that A and S denote analogous fine scale matrices as defined in Eq. @, 
except that averaged coefficients are used in the construction. To generate the offline space we 
then choose the smallest Mas eigenvalues from Eq. ([8]) and form the corresponding eigenvectors 
in the space of snapshots by setting ipf 1 = V . \I/^'0j nap (for k = 1, . . . , M off ), where \l/^ f are the 
coordinates of the vector \P£ ff • We then create the offline matrix 

R s =[<,..., <J 

to be used in the online space construction. 

For a given input parameter, we next construct the associated online coarse space VJ n (p) for 
each [i value on each coarse subdomain. In principle, we want this to be a small dimensional 
subspace of the offline space for computational efficiency. The online coarse space will be used 
within the finite element framework to solve the original global problem, where a continuous or 
discontinuous Galerkin coupling of the multiscale basis functions is used to compute the global 
solution. In particular, we seek a subspace of the offline space such that it can approximate any 
element of the offline space in an appropriate sense. We note that at the online stage, the bilinear 
forms are chosen to be parameter-dependent. Similar analysis (see [22]) motivates the following 
eigenvalue problem in the offline space: 

A°V)*r = \?£r(iM)9?, (9) 



where 

A° n (/i) = [a°» m „] = J K (x; /i) V^f • VC ff = i§fAO*)^ 

S°» = [s on (/iW] = J k(x; rilfflf = R^S(n)R om 

and k(x; //) and «(x; /i) are now parameter dependent. To generate the online space we then choose 
the smallest M on eigenvalues from Eq. (|9]) and form the corresponding eigenvectors in the offline 
space by setting ip™ = £\. ^™-ip° n (for k = 1, . . . , M on ), where \PS are the coordinates of the 
vector \I/£ n . 

5.2. Global coupling 

3.2.1. Continuous Galerkin coupling 

In this subsection we aim to create an appropriate solution space and variational formulation 
that is suitable for a continuous Galerkin approximation of Eq. (|5]). We begin with an initial coarse 
space V^ imt (/i) = spanjxj}^ (we use N v to denote the number of coarse vertices), where the x% 
are the standard multiscale partition of unity functions defined by 

- div (k(x; fj) Vxi) = K e uJi (10) 

Xi = 9i on dK, 

for all K E tUi, where g^ is assumed to be linear. Referring back to Eq. Q (for example), we note 
that the summed, pointwise energy k required for the eigenvalue problems will be defined as 



AT, 
K = K 

i=l 



E^ 2 i v ^i 2 - 



We then multiply the partition of unity functions by the eigenfunctions in the online space V£* to 
construct the resulting basis functions 

C = X^T' m for 1 < i < N v and 1 < k < M%, (11) 

where M^ denotes the number of online eigenvectors that are chosen for each coarse node i. We 



note that the construction in Eq. ( |TT| ) yields inherently continuous basis functions due to the mul- 
tiplication of online eigenvectors with the initial (continuous) partition of unity. This convention 
is not necessary for the discontinuous Galerkin global coupling, and is a focal point of contrast 
between the respective methods. However, with the continuous basis functions in place, we define 



the continuous Galerkin spectral multiscale space as 

V GG (fi) = span{^ CG : 1 < i < N v and 1 < k < M%}. (12) 

Using a single index notation, we may write V^ G (/i) = spanj^f }^, where N c denotes the total 
number of basis functions that are used in the coarse space construction. We also construct an 
operator matrix Rq = [ip CG , • • • , ip^*] (where ip CG are used to denote the nodal values of each 
basis function defined on the fine grid), for later use in this subsection. 

Before introducing the continuous Galerkin formulation, we recall that the parameter n is used 
to denote a solution that is computed at a previous iteration level (see Eq. ([5])). In turn, to update 
the solution at the current iteration level we seek m^(i; /i) = ^ Ciipf G (x; jjl) G V^ g such that 



CG/CG 



a [u 



ms i 



v,ti) = {f,v) for all veV GG , (13) 



where a (u, v; fi) — / k{x\ij)Vu- Vvdx, and (/, v) — \ fvdx. We note that variational form 

Jd J d 

in ( [13] ) yields the following linear algebraic system 

A U GG = F , (14) 

where U GG denotes the nodal values of the discrete CG solution, and 

Ao(ji) = [ajj] = [ k(x; ft) V^ G ■ V^ G dx and F = [fr] = [ M° dx. 
Jd Jd 

Using the operator matrix Rq, we may write A (fi) = RqA(^)R^ and F = R F, where A(/j,) 
and F are the standard, fine scale stiffness matrix and forcing vector corresponding to the form in 
Eq. ( fT3| ). We also note that the operator matrix may be analogously used in order to project coarse 
scale solutions onto the fine grid. 

3.2.2. Discontinuous Galerkin coupling 

One can also use the discontinuous Galerkin (DG) approach (see also flU [T51 EH) to couple 
multiscale basis functions. This may avoid the use of the partition of unity functions; however, 
a global formulation needs to be chosen carefully. We have been investigating the use of DG 
coupling and the detailed results will be presented elsewhere, see lfT8l . Here, we would like to 
briefly mention a general global coupling that can be used. The global formulation is given by 

a DG (n, v; /i) = f(v) for all v = {v K e V K }, (15) 

10 



where 

a DG (u,v;/i) = ^a£ G (w,v;/i) and f(v) = ^2 f v * dx i ( 16 ) 

K K ^ K 

for all u = {uk}, v = {v^} with K being the coarse element depicted in Figure [I] Each local 
bilinear form a^ G is given as a sum of three bilinear forms: 

a° G (w, v; fi) := a K (u, v; /i) + r K (u, v; //) + p K (u, v; //), (17) 

where a K is the bilinear form, 

o,k(u,v ; /i) := / Kk(^; /-OVwk • Vvj^dx, (18) 

where k k (x; //) is the restriction of k(x; //) in K; the r K is the symmetric bilinear form, 
r K (u, v,(jl):= ^ — / «jg(a?; fi) ( q-^( v k ~ v K >) + q~^( u ^' ~ m k) J ds, 

EcdK E J E V K K / 

where «e(x; fi) is the harmonic average of k(x; fi) along the edge E,Ie = 1 if E is on the boundary 
of the macrodomain, and Ie = 2 if E is an inner edge of the macrodomain. Here, K' and K are two 
coarse-grid elements sharing the common edge E; and px is the penalty bilinear form, 

ii r 

p K (u,v;n) := y^ -tt—^e / k e (x; fJ,)(u K > - u K )(v K < -v K )ds. (19) 

EcdK * * Jh 

Here h E is harmonic average of the length of the edge E and E', Se is a positive penalty parameter 
that needs to be selected and its choice affects the performance of GMsFEM. One can choose other 
eigenvalue problems within the DG framework. See [18J. 

As mentioned before that for Discontinuous Galerkin formulation, the support of basis func- 
tions are coarse element K as depicted in Figure [T] Besides, the inherent unconformal property of 
DG formulation determines the removal of the partition of unity functions while constructing basis 
functions in Equation ( fTTj ). Similarly, we can obtain the discontinuous Galerkin spectral multiscale 
space as 

K° G (/i) = span{^ DG : 1 < k < M*}, (20) 

For every coarse element K. 

Using the same process in the continuous Galerkin formulation, we can obtain an operator 
matrix constructed by the basis functions of V^ G (/i). For the consistency of the notation, we 



11 



denote the matrix as R , and Rq = Upx°, ■ ■ ■ , , 0° G ] • Recall that N c denote the total number of 
coarse basis functions. 

Solving the problem ([T]) in the coarse space V^ G (/i) using the DG formulation described in 
Equation ( fT5| ) is equivalent to seeking w^ G (x; fi) = ^V Cj^ 4 DG (x; /i) G V® G such that 



a UG «\ v- //) = f(v) for all v G V™, (21) 



where a DG («, v; /x) and /(t>) are defined in Equation p6] ). Similar as the CG case, we can obtain a 
coarse linear algebra system 

A U™ = F , (22) 

where Uq G denotes the discrete coarse DG solution, and 

Ao(ji) = R A(fi)RT and F = R F, 

where A{ji) and F are the standard, fine scale stiffness matrix and forcing vector corresponding to 
the form in Eq. ( fT6] ). After solving the coarse matrix, we can use the operator matrix R to obtain 



the fine-scale solution in the form of RqUq G . 



4. Numerical Results 



In this section we solve the nonlinear, elliptic model equation given in Eq. ([2]) using both the 
continuous (CG) and discontinuous Galerkin (DG) GMsFEM formulations described in Sect. |3j 
More specifically, we consider the equation 

-div(e K{x)u{x) Vu{x)) = f in D (23a) 

u = on dD, (23b) 

where the general coefficient from (Q is taken to be k(x; u) = e K ( x M x ). For the coefficient k(x), 
we consider the high-contrast permeability fields as illustrated Fig. [2] Fig. 2(a) represents a field 



whose high-permeability values are randomly assigned, while the field in Fig. 2(b) has a different 
channelized structure with fixed maximum values. We use a source term / = 0.1, and solve the 
problem on the unit two-dimensional domain D = [0, 1] x [0, 1]. 



To solve Eq. (23 ) we first linearize it by using a Picard iteration. In particular, for a given initial 
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Figure 2: High-contrast permeability fields 



guess u° we solve 



-div(e Kix)un{x) Vu n+1 (x)) = / in D 

ondD, 



u 



n+l 



(24a) 
(24b) 



for n > 0. 

In our simulations, we take the initial guess u° 



0, and terminate the iterative loop when 



\\A(u 



n+V 



II: 



n+l 



b || < 5 || b ||, where 5 is the tolerance for the iteration and we select 5 



10" 



We note that A and b correspond to the linear system resulting from either the CG or DG global 
formulations. In particular, we solve the problem as follows: 



A(u n )u 



n\ n+l 



b for n = 0, 1, 



(25) 



We note that since u n and u n+l will not necessarily be computed in coarse spaces of the same 
dimension, we cannot directly use the residual criterion listed above. Actually, we use the Galerkin 
projection of the fine solution to the corresponding coarse space to calculate the residual error from 
above. 



Remark 1. In this section we will consider two types of coefficients k(x) to be used in Eq. p3] ). 
We recall that throughout the paper we have used an auxiliary variable \x = u n to denote the 
solution dependence of the nonlinear problem. As such, we have referred to the model equation 
as parameter-dependent while describing the iterative solution procedure. Consequently, we are 
careful to introduce (and distinguish) a related case where we use a "physical " parameter //* for 
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Figure 3: Decomposition of permeability field 2(b) 



the purpose of constructing afield of the form k(x) = jj, p ki(x) + (1 — /j, p )k 2 (x). See Fig. \3jfor 
an illustration of K\{x) and K2{x). We note that the coefficient will be constructed by summing 
contributions that depend on the physical parameter pP, in addition to the auxiliary parameter 



dependence from the iterative form. In Subsect. 4.1 we use afield that does not depend on pP, and 



in Subsect. 4.2 we use afield that does depend on pP. 



4.1. Parameter-independent permeability field 

In the following simulations we first generate a snapshot space, use a spectral decomposition to 
obtain the offline space, and then for an initial guess apply a similar spectral decomposition to ob- 
tain the online space. We recall that in order to construct the snapshot space we choose a specified 
number of eigenfunctions (denoted by M snap ) on either a coarse neighborhood or coarse element 
depending on whether we use continuous (CG) or discontinuous Galerkin (DG) global coupling, 
respectively. In our simulations, we select the range of solutions [u 



mm; ^max 



that correspond to 

solving the fine scale equation using a source term that ranges from / G [0.1, 1]. For the first set 
of simulations we divide the domain [u min , u max ] into N s — 1 equally spaced subdomains to obtain 
N s discrete points U\, . . . , un s - For these simulations we fix a value of N s = 9. 



For either formulation, we solve a localized eigenvalue problems as defined in Subsect. 3.1 



for each point Uj on a coarse neighborhood and keep a specified number of eigenfunctions. For 
example, in the CG case we keep Z max = 3 snapshot eigenfunctions, and this construction leads 
to a local space of dimension M t 



snap 



I max xiV s = 3x9 = 27. In the DG case, we adaptively 
choose the number of eigenfunctions based on a consideration of the eigenvalue differences. In 
the offline space construction we fix u as the average of the previously defined fixed snapshot 
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values. We then solve the offline eigenvalue problem and construct the offline space by keeping 
the eigenvectors corresponding to a specified number of dominant eigenvalues. At the online stage 
we use the initial guess u° = in order to solve the respective eigenvalue problem required for the 
space construction. We note that the size of our online space and the associated solution accuracy 
will depend on the number of eigenvectors that we keep in the online space construction. 

In the CG formulation, we recall that the online eigenfunctions are multiplied by the corre- 
sponding partition of unity functions with support in the same neighborhood of the respective 



coarse node. We then solve Eq. (23 ) iteratively within the online space. In particular, for each it- 



eration we update the online space and solve the equation Eq. ( |23| ) using the previously computed 
solution. 

In the simulations using the CG formulation we discretize our domain into coarse elements of 
size H — 1/10, and fine elements of size h = 1/100. The results corresponding to the permeability 



fields from Figs. 2(a)| and [2(b)| are shown in Tables [T] and [2} respectively. The first column shows 



the dimension of the online solution space, and the second column shows the eigenvalue A* which 
corresponds to the first eigenfunction that is discarded from space enrichment. We note that this 
eigenvalue is an important consideration in error estimates of enriched multiscale spaces dl22~l0 . 
As a formal consideration, we mention that the error analysis typically yields estimates of the 
form \\u — u ms \\ ~ 0(H y X*) when the dominant eigenvalues are taken to be small. The next 
two columns correspond to the L 2 -weighted relative error \\u — ^msll^m / ||w|| L 2( £) ) xl00% and 
energy relative error \\u — Wms||#im) / IMI#im) x 100% between the GMsFEM solution u ms and 
the fine-scale solution u. We note that as the dimension of the online space increases (i.e., we keep 
more eigenfunctions in the space construction), the relative errors decrease accordingly. As an 
example, for the field in Fig. 2(a) we encounter L 2 relative errors that decrease from 1 .43 — 0.24%, 



and energy relative errors that decrease from 16.12 — 6.85% as the online space is systematically 
enriched. In the tables, analogous errors between the online GMsFEM solution and the offline 
solution are computed. The dimension of the offline space is taken to be the maximum dimension 
of the online space. We note that in this case the Picard iteration converges in 4 steps for all 
simulations. In Fig. [4] we also plot the fine and coarse-scale CG solutions that correspond to the 
field in Fig. 2(b)| We note that the fine solution, and the coarse solutions corresponding to the 



largest and smallest online spaces are nearly indistinguishable. 

We also illustrate the relation between the energy relative errors and A* in Fig. [5] for the same 
permeability fields considered above. From the plots in Fig. |5J we see that the energy relative error 
predictably decreases as A* decreases, thus following the appropriate error behavior. 

In order to solve the model problem using the DG formulation, we note that the space of snap- 
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dim(K C „ G ) 


A* 


GMsFEM Relative Error (%) 


Online-Offline Relative Error (%) 


Ll(D) 


mm 


Ll(D) 


H&D) 


319 


0.0021 


1.43 


16.12 


1.25 


16.33 


497 


0.0010 


0.69 


11.71 


0.48 


10.66 


770 


3.36 x 10" 4 


0.40 


9.13 


0.20 


7.30 


1043 


1.06 x 10" 4 


0.31 


7.76 


0.09 


4.43 


1270 


— 


0.24 


6.85 


0.00 


0.00 



Table 1: CG relative errors corresponding to the permeability field in Fig. 2(a) 



dim(K C n G ) 


A* 


GMsFEM Relative Error (%) 


Online-Offline Relative Error (%) 


Ll(D) 


HkiP) 


Ll(D) 


Hk(D) 


316 


0.0026 


1.36 


15.28 


1.18 


15.74 


482 


0.0010 


0.71 


11.89 


0.51 


11.17 


722 


3.18 x 10~ 4 


0.43 


9.53 


0.22 


7.77 


996 


1.02 x 10" 4 


0.33 


8.02 


0.11 


4.72 


1236 


— 


0.26 


7.05 


0.00 


0.00 



Table 2: CG relative errors corresponding to the permeability field in Fig. 2(b) 



shots is constructed in a slightly different fashion. In this case, the selection of eigenvectors hinges 
on a comparison between the difference of consecutive eigenvalues resulting from the localized 
computations. In contrast to the CG case, the initial number of eigenf unctions (call this number 
Z^j t ) used in the snapshot space construction are adaptively chosen based on the relative size of 
consecutive eigenvalues. For the results corresponding to the DG formulation, we note that two 
configurations for the snapshot space construction are used. In particular, we consider a case when 
the original number of eigenf uctions l? it are used in the construction, and a case when l^ ax = 1^+3 
are used in the construction. 



CG Fine Solution 




CG GMsFEM 



= 31/ US 



0.2 0.4 0.6 0.8 1 




CG GMsFEM: 



= 1270 xio 



0.2 0.4 0.6 0.8 




0.2 0.4 0.6 0.8 



Figure 4: Comparison of fine and coarse CG solutions correpsonding to Fig. 2(b) 
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(a) Corresponds to Fig. 2(a) 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 
A* xlCT 3 



(b) Corresponds to Fig. 2(b) 



Figure 5: Relation between the first discarded eigenvalue and the CG relative energy error; permeability from Fig. 2(a) 



(left), permeability from Fig. 2(b) (right) 



In the simulations using the DG formulation, we partition the original domain using a coarse 
mesh of size H = 1/10, and use a fine mesh composed of uniform triangular elements of mesh 
size h = 1/100. The numerical results for permeability fields 2(a) and 2(b) are represented in 



Tables [3] and HI respectively. The first column shows the dimension of the online space, the second 
column represents the corresponding eigenvalue(A*) of the first eigenfunction discarded from the 
online space, and the next two columns illustrate the interior energy relative error (E int ) and the 
boundary energy relative error (Eg) between the fine scale solution and DG GMsFEM solution. 
The errors between the offline and online solutions are offered in the final two columns. We 
note that as the dimension of the online space increases (i.e., we keep more eigenfunctions in 
the space construction), the relative errors decrease accordingly. For example, the DG solution 



corresponding to Fig. 2(a) yields interior relative energy errors that decrease from 55.08 — 34.86%, 
and boundary relative energy errors that decrease from 8.94 — 6.40%. We note that in this case the 
Picard iteration converges in 4 or 5 steps for all simulations. In Fig. [6] we also plot the fine and 



coarse DG solutions that correspond to the field in Fig. 2(b) We note that the fine solution and the 
coarse solution corresponding to the smallest online space show some slight differences. However, 
the discrepancies noticeably diminish when the coarse DG solution is computed within the largest 
online space. As in the CG case, we also illustrate the relation between the DG interior errors and 
A* in Fig.[7J From the plots in Fig. [7} we see that the relative errors decrease as A* decreases, again 
following the expected error behavior. 

Remark 2. When solving the nonlinear equation using the discontinuous Galerkin approach, we 
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dim(yD°) 


A* 


GMsFEM Relative Error (%) 


Online-Offline Relative Error (%) 


Eint 


Eg 


-E-int 


Eg 


271 


1.53 x 10" 4 


55.08 


8.94 


44.38 


8.43 


331 


1.24 x 10" 4 


36.59 


6.63 


10.05 


3.08 


466 


3.03 x 10" 5 


35.57 


6.56 


7.00 


1.67 


624 


1.72 x 10" 5 


34.90 


6.48 


2.12 


0.40 


716 


— 


34.86 


6.40 


0.00 


0.00 



Table 3: DG relative errors corresponding to the permeability field in Fig. 2(a) snapshot space uses l^ h eigenfunctions 



dim(K° G ) 


A* 


GMsFEM Relative Error (%) 


Online-Offline Relative Error (%) 


^int 


Eg 


Ei„t 


Eg 


270 


1.56 x 10~ 4 


56.29 


10.30 


46.37 


9.75 


331 


1.05 x 10" 4 


36.72 


6.71 


9.54 


3.32 


444 


3.12 x 10" 5 


35.67 


6.56 


6.48 


1.67 


582 


1.21 x 10~ 5 


35.06 


6.48 


2.14 


0.41 


663 


— 


35.03 


6.48 


0.00 


0.00 



Table 4: DG relative errors corresponding to the permeability field in Fig. 2(b) , snapshot space uses l^ it eigenfunctions 



use different penalty parameters for fine-grid problem and coarse-grid problem (refer back to 



Subsect. 3.2.2). However, we observe that for different coarse penalty parameters that yield a 



convergent solution, the number of iterations and the relative errors (both interior and boundary) 
stay the same. 

Remark 3. Recall that we use the Galerkin projection of the previous coarse solution onto the 
current online space as the approximation of the previous coarse solution to obtain the terminal 
condition. If the coarse penalty parameter is changed, we should use the current coarse penalty 
parameter to construct the Galerkin projection. 



DG Fine Solution 




DG GMsFEM: dim! 



DG GMsFEM: dim(K ») = 663 



0.2 0.4 0.6 0.8 




Figure 6: Comparison of fine and coarse DG solutions correpsonding to Fig. 2(b) 
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g? 60 




(a) Corresponds to Fig. 2(a) 



(b) Corresponds to Fig. 2(b) 



Figure 7: Relation between the first discarded eigenvalue and the DG relative interior energy error; permeability from 



Fig. 2(a) (left), permeability from Fig. 2(b) (right) 2(a) and 2(b) 



dim(K° G ) 


A* 


GMsFEM Relative Error (%) 


Online-Offline Relative Error (%) 


^int 


E a 


^int 


E d 


381 


1.47 x 10~ 4 


37.34 


7.42 


22.80 


6.00 


440 


1.54 x 10~ 4 


35.92 


6.16 


20.07 


4.36 


707 


9.54 x 10~ 5 


32.80 


5.29 


13.64 


2.90 


958 


2.71 x 10~ 5 


29.44 


5.48 


4.80 


0.98 


1352 


— 


28.98 


5.39 


0.00 


0.00 



Table 5: DG relative errors corresponding to the permeability field in Fig 2(b), snapshot space uses l^. dx = l{ 
eigenfunctions 



K 

init 



We observe from Tables [T]j4j that the offline spaces for DG formulation are much smaller than 
those obtained through CG formulation. As a result, in Table [5] we use more eigenfunctions (more 
specifically, we set Z^ ax = l^ it + 3) in the snapshot space construction to yield a larger offline space. 



For these examples, we use the permeability field from Fig. 2(b) Due to the increase of the offline 
(and corresponding online) space dimensions, we see more accurate results than those offered in 
TableS 

4.2. Parameter-dependent permeability field 

For the next set of numerical results, we consider solving the nonlinear elliptic problem in 
Eq. ( |23| ) with a coefficient of the form k(x, u, /i p ) = exp [(fJ'n^x) + (1 — /i p )K 2 (x)) u(x)}. For 



K\{x) and k 2 (x) we use the fields shown in Fig. 3(a) and 3(b), respectively. As for the parameter 



dependent simulation, we are careful to distinguish the difference between the auxiliary parameter 
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dim(K C „ G ) 


A* 


GMsFEM Relative Error (%) 


Online-Offline Relative Error (%) 


Ll(D) 


H^D) 


Ll(D) 


H&D) 


309 


0.0027 


1.30 


14.89 


1.10 


15.32 


492 


0.0010 


0.59 


10.82 


0.39 


9.76 


580 


6.76 x 10" 4 


0.45 


9.55 


0.24 


7.92 


728 


3.33 x 10" 4 


0.34 


7.87 


0.12 


5.23 


991 


— 


0.28 


6.74 


0.00 


0.00 



Table 6: CG relative errors corresponding to the parameter-dependent field constructed from Fig. 3(a) and 3(b) 



dim(K° G ) 


A* 


GMsFEM Relative Error (%) 


Online-Offline Relative Error (%) 


^int 


E a 


Emt 


E d 


300 


1.02 x 10~ 4 


37.56 


7.94 


10.15 


3.16 


313 


6.25 x 10~ 5 


37.55 


7.81 


10.00 


2.85 


403 


2.58 x 10~ 5 


36.81 


7.35 


5.83 


1.38 


497 


1.22 x 10~ 5 


36.37 


7.21 


0.84 


0.10 


517 


— 


36.36 


7.21 


0.00 


0.00 



Table 7: DG relative errors corresponding to the parameter-dependent field constructed from Fig. 3(a) and 3(b) 
snapshot space uses l(£ it eigenfunctions 



fj, = u n which is used to denote a previous solution iterate, and a "physical" parameter jjP that is 
used in the construction of a new permeability field. We take the range of fjP to be [0, 1], and use 
three equally spaced values in order to construct the snapshot space in this case. We use the same 
interval from the previous results, yet use four equally spaced values in this case. In 



u 



mm j Mmax 

particular, we use the pairs (uj, /J}), where 1 < j '< 4, and 1 < I < 3 as the fixed parameter values 
for the snapshot space construction. At the online stage we use the initial guess u° = and a fixed 
value of /j? = 0.2 while solving the respective eigenvalue problem required for the continuous or 
discontinuous Galerkin online space construction. 

In Table[6]we offer results corresponding to the CG formulation, and in Tables [7] and [8] we offer 
results corresponding to the DG formulation. In all cases we encounter very similar error behavior 
compared to the examples offered earlier in the section. In particular, an increase of the dimension 
of the online space yields predictably smaller errors, and smaller values of A* correspond to the 
error decrease. And while it suffices to refer back to related discussions earlier in the section, we 
emphasize that this distinct set of results serves to further illustrate the robustness of the proposed 
method. In particular, we show that the solution procedure allows for a suitable treatment of 
nonlinear problems that involve auxiliary and physical parameters. 
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dim(K? G ) 


A* 


GMsFEM Relative Error (%) 


Online-Offline Relative Error (%) 


Eint 


Eg 


Eint 


Eq 


300 


2.13 x 10" 4 


37.59 


7.94 


22.54 


6.40 


440 


1.54 x 10" 5 


35.78 


5.92 


18.89 


3.74 


668 


7.69 x 10" 5 


32.54 


5.39 


11.62 


2.58 


902 


1.51 x 10" 5 


30.23 


5.29 


3.87 


1.06 


1093 


— 


29.88 


5.29 


0.00 


0.00 



Table 8: DG relative errors corresponding to the parameter-dependent field constructed from Fig. 3(a) and 3(b) 
snapshot space uses I 



K 
max 



^init + ^ eigenfunctions 



5. Concluding Remarks 

In this paper we use the Generalized Multiscale Finite Element (GMsFEM) framework in order 
to solve nonlinear elliptic equations with high-contrast coefficients. In order to solve this type of 
problem we linearize the equation such that upscaled quantities of previous solution iterates may 
be regarded as auxiliary coefficient parameters in the problem formulation. As a result, we are able 
to construct a respective set of coarse basis functions using an offline-online procedure in which 
the precomputed offline space allows for the efficient computation of a smaller-dimensional online 
space for any parameter value at each iteration. In this paper, the coarse space construction involves 
solving a set of localized eigenvalue problems that are tailored to either continuous Galerkin (CG) 
or discontinuous Galerkin (DG) global coupling mechanisms. In particular, the respective coarse 
spaces are formed by keeping a set of eigenfunctions that correspond to the localized eigenvalue 
behavior. Using either formulation, we show that the process of systematically enriching the coarse 
solution spaces yields a predictable error decline between the fine and coarse-grid solutions. As a 
result, the proposed methodology is shown to be an effective and flexible approach for solving the 
nonlinear, high-contrast elliptic equation that we consider in this paper. 
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